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ABSTRACT 

The orbital distributions of currently observed extrasolar giant planets allow marginally stable orbits 
for hypothetical, terrestrial planets. In this paper, we propose that many of these systems may not have 
additional planets on these "stable" orbits, since past dynamical instability among giant planets could 
have removed them. We numerically investigate the effects of early evolution of multiple giant planets 
on the orbital stability of the inner, sub-Neptune-like planets which are modeled as test particles, 
and determine their dynamically unstable region. Previous studies have shown that the majority of 
such test particles are ejected out of the system as a result of close encounters with giant planets. 
Here, we show that secular interactions can remove test particles far from giant planets as well, and 
that such test particles tend to merge with the central stars. Our results indicate that, unless the 
dynamical instability among giant planets is either absent or quiet like planet-planet collisions, most 
test particles within the orbits of giant planets at a few AU may be gone. We find a good agreement 
between our numerical results and the secular theory, and present a semi-analytical formula which 
estimates the dynamically unstable region of the test particles just from the evolution of giant planets. 
Finally, our numerical results agree well with the observations, and also predict the existence of hot 
rocky planets in eccentric giant planet systems. 
Subject headings: planetary systems; methods: N-body simulations 



1. INTRODUCTION 

The Solar System has a rich planetary architecture 
with four inner terrestrial planets and four outer giant 
planets. Whether such a planetary system is common 
or not has been a long-standing question in astronomy. 
The recent observations started revealing that multiple- 
planet systems are common, espe cially among close- 
in sub-Neptune-like planets (e.g., iMavor et al.l 120111 : 
iBatalha et al.ll2012t ). However, we are yet to find sys- 
tems similar to ours. Future observations would be able 
to tell us whether some of these systems with multiple 
sub-Neptune-like planets are accompanied by giant plan- 
ets on further orbits, and/or known giant-planet systems 
have terrestrial planets interior to their orbits. 

In the meantime, we can explore whether the currently 
observed planetary systems allow the existence of ad- 
ditional planets. Such an attempt has been made by 
investigating the orbital stability of hypothetical terres- 
trial planets in the observed giant-planet systems (e.g . 
Gehman et al.lll996t| Jones et al.ll2001k iNqble et , al.l l2002l: 



Menou fe Tabachnikl 120031: lAsghari et all I2004D These 
studies showed that some of the observed systems could 
indeed host dynamicall y stable Earth-like planets in their 
habitable zones (HZs). Menou fc Tabachnikl (|2003f ) mod- 
eled terrestrial planets as test particles, and investigated 
dynamical habitability of these planets in 85 then-known 
exoplanetary systems. They found that about 25% of 
their samples could host habitable, Earth-like planets. 
These studies, however, do not take account of the ef- 
fects of potential early evolution of giant-planet systems 
on formation and/or survival of terrestrial planets. 
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The migration and dynamical instability of giant plan- 
ets are two major threats to formation an d stability 
of Ear th-like planets interior to their orbits. lArmitaga 
(|2003h pointed out the potentially hazardous effect of 
giant planet migration on formation of terrestrials, and 
suggested that these planets preferentially exist in sys- 
tems where massive giants did not migrate significantly. 
Various numerical simulations, however, are generally 
agreed that terrestrial planet formation is not necessar- 
ily prevented by giant planet migration, when eccen- 
tricity excitation timescales for (proto-)terrestrial plan- 
ets are long c ompared to migration timescales of giant 
planets (e.g.. iMandell fc Sigurdssonl 120031 : iLufkin et al.l 
120061 : iRavmond et al.ll2006f ). These less massive planets 
can survive by either being shepherded interior to the 
migratin g giant planet, or sca ttered exterior to its or- 
bit (e.g.. iFogg fc Nelson! 120051 ) . However, the survival 
rate of less massive planets strongly depe nds on the mi- 
grati o n speed of the mas s ive one (e.g., iWard fc Hahnl 
1991 iTanaka fc Idal fl999l: lEdgar fc Artvmowiczl \200W. 



Mandell fc Sigurdssonl (12001 " showed that the survival 
percentage of planets decreases from 40% to 15% when 
the migration timescale for a giant planet increases from 
0.5 to 2Myr. 

How the dynamical instability of giant planets af- 
fects terre strial planet formation was also investi- 
gated (e.g . , Theb ault et al.l 2002 ; iLevison fc Agnorll200l: 
IVeras fe Armitagd I2005L 120061 )! iVeras fc Armitagd 
( 200l ~ showed that, when a moderately eccentric giant 
planet is left at a few AU, over 95% of test particles are 
ejected — i.e., terrestrial planet formation is more dif- 
ficult in such a system. However, their study focused 
on the stability of planets in the HZs (0.75 - 1.25 AU). 
Although their follow-up study has extended the orbital 
radii down to 0.5 AU. the recent observations have found 
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many sub-Neptune-sized planets within the orbital ra- 
dius of Mercury (about 0.4 AU). Thus, we must take ac- 
count of much closer-in planets to discuss the effects of 
dynamical instability among giant planets on the survival 
of these short-period, low-mass planets. 

In this paper, we study the orbital evolution of ter- 
restrial planets when giant planets become dynamically 
unstable, and show that terrestrial planets far away from 
giants can also be removed. The numerical methods and 
initial conditions of our simulations are described in Sec- 
tion 2. In all of our simulations, we approximate terres- 
trial planets as test particles, and follow their evolution 
as three Jupiter-mass planets become dynamically un- 
stable. As we show in Section 3, many of these planets 
merge with the central stars, rather than being ejected 
out of the systems. In Section 13.11 we compare the nu- 
merical results with the secular theory, and show that 
mergers are largely due to secular interactions. In Sec- 
tion !3.21 we present a convenient formula which estimates 
the dynamically unstable region of terrestrial planets 
from the orbital evolution of a giant-planet system alone. 
In Section 13.31 we discuss how the fates of test particles 
depend on the dynamical outcomes of giant planets. In 
Section [3.4[ we investigate the importance of past evolu- 
tion of giant planets on the orbital stability of test par- 
ticles. In Section 13.51 we compare our results with the 
distributions of observed exoplanetary systems. Finally, 
we address the directions of future studies, and summa- 
rize our work in Section 4. 

2. METHOD 

We perform numerical simulations of planetary sys- 
tems by using a p ublicly available N-body code Mercury 
(|Chambersl 119991) . For all the simulations presented in 
this paper, we use the Bulirsch-Stoer method with a tol- 
erance parameter of 10~ 12 to ensure the high accuracy. 

The standard planet formation theory predicts that gas 
giant planets are formed beyond the ice line. However, 
there is no consensus about how many such planets being 
formed per a typical protostellar disk. Since we are inter- 
ested in the effects of dynamical instability among giant 
planets on survival or fo rmation of terrestr i al pla nets, 
we follow the approach of iVeras fc Armitaed (|2005l ) and 
use three giant-planet systems. This is also the smallest 
number of plan ets for which there is no analytical sta - 
bility limit (e.g.. lMarchal fe Bozid[l98llGladmanlll993f ). 
Also, we implicitly assume that giant planets did not mi- 
grate significantly in the disk. Such simulations would 
not be able to recover the orbital distribution of ob- 
served exoplanets, especially that within ~ 1 AU (e.g., 
lAdams fc Laughlinl |200H iChatteriee et all 120081) . The 
closer-in plan ets could be formed either via planet mi- 
gration ("e.g.. iLin et al.lll996t lAdams fc Laughlinl 120031 : 
iMoorhead fc A dams 2008; Matsumura et al.ll2009f l or via 
dynamical evolution of giants co mbined with tidal inter- 
actions with the cent r al star (e.g. iFabrvckv fc Tremaine 
20071 IWuet all 120071: INagasawa et all 120081; iNaoz et al 
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IWu fc Lithwickil201in . We ignore these effects here, 



and focus on the importance of the dynamical instability 
among giant planets on the inner planetary systems. 

We assume that all three giant planets have a Jupiter 
mass, and the initial semimajor axes of 3.8, 5, and 6.5 
AU, respectively. Such an orbital distribution generally 
leads to a dynamical instability within a short period 



of time. Naturally, the actual planet formation is un- 
likely to lead to such quasi-unstable systems. However, 
once the instability occurs, the final outcome of dynam- 
ical evolution of planets will not be largely affected by 
the initial distribution of planets, as long as the total en- 
ergy and the angular momentum do not differ m uch (e.g., 
Uuric fc Tremainel 120081 : INagasawa et alj 12008ft . There- 
fore, our choice of the orbital distribution of planets is 
for the computational convenience. We further assume 
equal-mass giant planets, beca use they lead to the small- 
est mean pericenter distance (jVeras fc Armitagel I2006T ) 
and thus give the most severe condition for the survival 
of terrestrial planets. Note that the instability among 
equal-mass planets tends to result in the high eccentric- 
ities of sur vived planets compa red to unequal-mass sys- 
tems (e.g.. lFord fc Rasioll2008ft . Thus, we do not intend 
to match the eccentricity distribution of giant planets 
with the observed one either. However, since we are only 
interested in the range of the stable region, and not the 
rate of stable test particles, this assumption would not 
change our conclusion. We also assume that initial ec- 
centricities of Jupiters are zero, and they have non-zero, 
but small initial inclinations of 0.001, 0.005, and 0.003 
degrees, respectively. 

Furthermore, we model terrestrial planets as test parti- 
cles in this work, which is a reasonable assumption since 
the mass ratio of Earth to Jupiter is small (Me/Mj ~ 
3 x 10" 3 ). We have 11 test particles over 0.1 to 3 AU. 
All of them initially have circular and coplanar orbits. 

We perform 40 scattering simulations, where we only 
change the initial orientations of giant planets' orbits. 
Each test particle has a randomly assigned orbital phase, 
but their distribution is common for all the simulations. 

3. RESULTS 

In Section 13.11 we present the overall results of our 
simulations. We find a good agreement between our re- 
sults and the expectations from the secular theory. In 
Section [3. 2\ we propose a formula which can predict the 
dynamically unstable region of test particles. In Sec- 
tion 13.31 we discuss the dependence of fates of test par- 
ticles on dynamical outcomes of giant planets. We show 
that the dynamical instability of giants would change the 
stable region of test particles in Section 13.41 Our results 
are compared with the observations in Section 13.51 

3.1. Fates of Test Particles depending on Minimum 
Orbital Distance of Jupiters 

We evolve each of 40 systems up to 10 Myr. When sys- 
tems have two Jupiters and more than one survived test 
particles at the end of the simulations, we run them up 
to 100 Myr. However, the overall results do not change 
significantly. 

An example of the orbital evolution of our model is 
shown in Fig. [TJ Three thick curves and eleven thin 
curves correspond to semi-major axes of Jupiter-mass 
planets and test particles, respectively. Dashed curves 
represent each particle's peri- and apo-center. This ex- 
ample shows a relatively "quiet" dynamical evolution of 
giant planets, where their orbital radii and eccentrici- 
ties do not change much — the outermost Jupiter-mass 
planet crosses orbits with inner two giants, and collides 
with the innermost Jupiter-mass planet around 2000 yr. 
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Fig. 1. — Orbital evolution of three Jupiter-mass planets (thick 
curves) and eleven test particles (thin curves). Note that one of 
the test particles hits the star before 100 yr and thus does not 
appear in this plot. Dashed curves indicate peri- and apo-centers 
of each particle. The orange, green, and blue colors of test particles 
indicate their fates of ejection, merger with the central star, and 
survival for 10 Myr, respectively. 




Fig. 3. — Same as Fig. [T] but without peri- and apo-centers, 
and with forced eccentricity contours which are estimated by the 
lowest-order secular theory with the orbital parameters of Jupiter- 
mass planets at each time. Orbits of merged test particles overlap 
with the high forced-eccentricity regions. 
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Fig. 2. — Orbital evolution of a merged and an ejected test par- 
ticle in Fig. [T] The 8th and 10th test particles from the star are 
shown with green and orange curves, respectively. Top, middle, 
and bottom panels show evolution of semimajor axes, eccentrici- 
ties, and inclinations, respectively. For the ejected test particle, 
both eccentricity and semimajor axis change within a short period 
of time due to strong gravitational interactions with Jupiter-mass 
planets. On the other hand, for the merged test particle, its semi- 
major axis stays constant while the eccentricity increases due to 
secular interactions, which leads to a small pericenter distance and 
an eventual merger with the central star. 

On the other hand, test particles appear to take vari- 
ous dynamical outcomes, largely depending on their dis- 
tances from Jupiters. The outermost test particle hits 
the central star immediately after the start of the sim- 
ulation (~ 36 yr), and thus is not shown in this plot. 
For the rest of the test particles, the outer test particle 
(orange curve) experiences orbital crossings with giants 
and is ejected out of the system, while the inner test par- 
ticles (blue curves) stay dynamically stable for 10 Myr. 
Over the intermediate radii, test particles which do not 
have orbital crossings with giant planets merge with the 
central star (green curves). 

When the test particles have orbital crossings with gi- 
ant planets, they tend to experience close encounters and 
strong interactions, which could dramatically change the 
orbits of test particles. On the other hand, even when the 
test particles do not have orbital crossings with giants, 



their orbits could still be perturbed by secular interac- 
tions. In Fig. [2j we highlight these differences between 
ejected and merged test particles. Since the ejected test 
particle (orange curves) has orbital crossings with giants, 
the semi-major axis changes rapidly while the eccentric- 
ity and inclination increase. On the other hand, for a 
merged test particle shown here (green curves), secular 
interactions dominate its evolution, and thus semi-major 
axis stays constant while the eccentricity increases. Fig- 
ures[T]and[2]indicate that, since secular interactions don't 
change a semi-major axis, the increase of an eccentricity 
leads to a smaller and smaller pericenter of a test particle 
and an eventual merger of the particle with the central 
star, rather than the ejection out of the system. 

To confirm the effects of secular interactions on or- 
bital evolution of test particles, we compare the expecta- 
tions from the secular theory with the simulated results. 
We adopt the lowest-order secular theory, and thus as- 
sume small eccentricities and inclinations for the planets. 
This is a rather crude assumption, because eccentricities 
and inclinations of giant planets can become very large 
as a result of the dynamical instability. Furthermore, 
since test particles are distributed over 0.1 — 3 AU in our 
model, their eccentricities must increase above 0.95 to hit 
the Sun-like central star. Nevertheless, as shown below, 
we find that the lowest-order secular theory predicts the 
fates of test particles fairly well. 

For the lowest-order secular theory, the evolutions of 
eccentricity and inclination are decoupled. Here, we fo- 
cus on the eccentricity evolution, because that deter- 
mines whether a test particle merges with the star or 
not. The secular solutions for eccentricity evolution con- 
sist of free and forced components. 

e sin w = e frco sin (At + /3) + e forccd sin n7 forcod (1) 
e cos w = e frco cos (At + j3) + ef orcc d cos n7 forccd (2) 

Here, e is the eccentricity, w is the longitude of pericen- 
ter, and (3 is the phase. The free component is character- 
ized by a free frequency A of the test particle's apsidal 
precession due to gravitational interactions with giant 
planets, while the forced component is characterized by 
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Fig. 4. — Fates of test particles depending on the initial semi- 
major axes of test particles that are normalized by the minimum 
pericenter distances (top panel) and the minimum semimajor axes 
of Jupiter-mass planets (bottom panel). The number of the ejected 
test particles increase toward and beyond a-tp.in/ &j min = b while 
most merged particles are within this limit. 
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Here, Vi = 'EjL-^Ajeji, where e,j is a component of an un- 
sealed eigen vector which corresponds to gi, and Aj de- 
pends on masses of a giant planet and the central star as 
well as the orbital r adii of the giant planet and a test par- 
ticle (for details, see lMurrav fc Dermottl[l999l ) . A similar 
equation holds for the cosine term. These equations in- 
dicate that the secular interactions become the strongest 
near resonance locations, where a free frequency of a test 
particle becomes comparable to one of the eigen frequen- 
cies of these giants A ~ gi. Since the efficient angular 
momentum exchanges lead to high forced eccentricities 
near secular resonances, we expect that test particles in 
these regions are prone to mergers. 

In Fig. [3l we compare the evolution of forced eccentric- 
ity contours with the orbital evolution shown in Fig. [TJ 
The orbits of merged test particles coincide with high 
forced eccentricity regions, while those of survived test 
particles are in small forced eccentricity regions. We have 
performed similar comparisons for all of our simulations, 
and found that the forced eccentricities of test particles 
are good indicators of removed and survived planets. 

If the secular interaction is a dominant cause of merged 
particles while the close encounters is a dominant cause 
of ejected particles, we would expect that the fates of 
test particles are divided clearly depending on whether 
test particles cross orbits with giant planets or not. We 
test this hypothesis in Fig. @] where the test particles' 
initial semimajor axes at p ^ n are compared with the min- 
imum pericenter qj, m in (top panel) and semimajor axes 
a>j,min (bottom panel) of giant planets. Note that, the 
"minimum" pericenter (or semimajor axis) of a Jupiter 
can be different for each test particle in the same sim- 
ulation, because this is chosen during each test parti- 
cle's survival time (i.e., if the test particle survived until 
the end of the simulation, the minimum semi-major axis 
would be that throughout the simulation). Surprisingly, 
we find that about half of the merged test particles have 
atp,in/qj,min > 1, and thus are potential candidates for 
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Fig. 5. — Comparisons of the mutual inclinations (top panels), 
angles between pericenters (middle panels), and minimum orbital 
distances (bottom panels). Left panels show the results when the 
pericenter of the innermost giant planet takes the smallest value 
Q J mini while right panels show the corresponding results at the 
end of the test particles' lifetimes (or at the end of the simulations 
for survived ones). The color codes are the same as Fig. [4] except 
that the merged test particles with at p m/lJ,min _ 1 are plotted 
in black instead of dark green. The middle panels show the cu- 
mulative histograms of the angles between pericenters A8 for each 
case corresponding to the top panels. The bottom panels show the 
minimum distance between orbits of test particles and the inner- 
most giant planets. This distance peaks near 10 Hill radii of the 
innermost giant planets for the survived test particles, and near a 
few Hill radii for both merged and crossed-then-merged particles. 

orbital crossings (see top panel). By checking the or- 
bital evolution of planets for each simulation, we find 
that some merged test particles have clear orbital cross- 
ings with giant planets, just like ejected test particles. 
However, most merged test particles have nearly con- 
stant semimajor axes until the merger times, and thus 
are expected to have negligible close encounters with gi- 
ant planets. The expectation that atp,m / Qj,min •_ 1 does 
not necessarily mean the close encounters, or even the 
orbital crossings, is confirmed by the bottom panel of 
Fig.|4l where atp. in I 'o-j,min > 1 indicates the obvious or- 
bital crossings. In this plot, few merged test particles 
exist beyond (Hp,inl 'a>J,min = F 

Fig. [5] shows that the "crossed-then-merged" (merged 
with a tp ^n/qj,min > 1) test particles actually avoid close 
encounters with giants due to large mutual inclinations 
between orbits, and/or a tendency toward an orbital 
alignment. From the top left panel, we find that, when 



is reached, the mutual inclination between orbits 



peaks near 0° for all cases except for the crossed-then- 
merged ones. The clear tendency toward large mutual 
inclinations for the crossed-then-merged cases become 
stronger at the end of the simulations (top right panel). 

For the middle panels, we introduce an angle between 
pericenters of two orbits AO. We use this angle as an 
indicator of the orbital alignment instead of the angles 
between argument of pericenters Aw, because Aw ~ 0° 
does not necessarily mean an alignment of pericenters 
depending on the mutual inclination. In both left and 
right panels, we find that survived test particles tend to 
be more aligned (A6 ~ 0°) than the others. The crossed- 
then-merged cases do not show such a strong tendency 
toward the alignment, but they avoid the misalignment 
(A6* ~ 180°) which would place the apocenter of the in- 
ner test particle close to the pericenter of the outer giant 
planet, and therefore could result in strong interactions 



between them. This tendency away from the misalign- 
ment becomes stronger toward the end of the simula- 
tions (middle right panel) , where the AO peaks near 70° . 
Therefore, these two sets of panels indicate that crossed- 
then-merged particles could have avoided close encoun- 
ters with giants due to large mutual inclinations and/or 
relatively well-aligned orbits. 

To check this point further, we compare the minimum 
distances between the orbits of a test particle and the 
innermost giant planet in the system (bottom panels). 
As expected, this distance is the largest for the survived 
test particles, and smaller for merged particles. For the 
ejected test particles, the distribution is broader, which 
indicates that many of these test particles started cross- 
ing orbits with giant planets before Qj^rnin 

was reached. 

For the survived test particles, the minimum orbital dis- 
tance peaks near 10 Hill radii of the giant planet, while 
for the merged and crossed-then-merged particles, this 
peaks near a few Hill radii. Although some crossed-then- 
merged test particles have very small minimum orbital 
distances (Ar < 1 Hill radius), the overall similarity 
of the distributions of merged and crossed-then-merged 
cases indicates that many crossed-then-merged particles 
have avoided close encounters with giant planets. 

From these results, we conclude that a tPt i n / 'o,j,min 
is the better indication of close encounters compared 
to dtp, in ll J, mint since the former shows the transition 
from merger-dominated to ejection-dominated cases near 
atp,in = a j, mim while the latter does not. The bottom 
panel of Fig. |4] confirms our expectation, and shows that 
the fates of test particles depend on distances from gi- 
ants — from the innermost region to outward, survived, 
merged, and ejected test particles are distributed. An- 
other notable point is a large fraction of the merged par- 
ticles compared to the ejected ones. Even in the "orbit 
crossing" cases with q^ mi „ < a tp j n (< aj >m in)i merger is 
domi nated over ejection . This is a very different outcome 
from lVeras fe A rmitagc (2005|), where most test particles 
are ejected. The difference is likely caused by the fact 
that IVeras fc~A rmitagc. (|2005h focused on the stability 
in HZ, while we are simulating test particles over a wider 
range of orbital radii. 

3.2. Estimating Dynamically Unstable Region of Test 
Particles 

In the last subsection, we have found that the secular 
interaction is a dominant cause of the mergers of test 
particles with the central stars. In this subsection, we 
test this hypothesis and propose a convenient formula 
which allows us to estimate the dynamically unstable re- 
gion of the test particles from the orbital evolution of 
giant planets alone. 

As we have seen in the last subsection, test parti- 
cles need to achieve e ~ 1 to hit the central stars. 
On the other hand, in the secular theory, the eccen- 
tricity is determined by the vector sum of the free and 
forced eccentricities. In the plane of (ecostu, esinnj), 
we can think that the free eccentricity vector completes 
a circular motion around the forced eccentricity vec- 
tor (e forcod costuf orced , ef orcod sinTZ7 forccd ) at the free fre- 
quenc y A while the forced vec tor itself is also moving 
(e.g., iMurrav fe Dermottl [l999h . Thus, we can expect 
that when the forced eccentricity of a test particle stays 
above 0.5 for the time longer than a secular timescale, 
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Fig. 6. — Comparison of the ratios of the total durations of 
the forced eccentricities of test particles being above 0.5 and the 
secular timescales. The forced eccentricities are estimated from 
the orbital parameters of Jupiters. This time ratio gives a good 
prediction of whether test particles being dynamically unstable or 
not. Most ejected and merged test particles exist beyond the ratio 
1, while most survived test particles are below this limit. A peak 
at 10 — 5 for survived test particles represents all the particles with 
Ai(e > 0.5)/t scc < 10~ 4 that do not achieve the estimated forced 
eccentricities above 0.5 throughout the simulations. 

the eccentricity may reach ~ 1. 

In Fig. [6l we test this hypothesis by calculating the 
forced eccentricity of a test particle from the orbital pa- 
rameters of giant planets and comparing the total dura- 
tion of ef orccd > 0.5 with a secular time: 

r S ec = min (2vr/A 2n/g. l ) . (4) 

Here, gi is an eigenfrequency of giants, and we choose 
the frequency that corresponds to the largest coeffi- 
cient Vi/{A — gi) in Eq. [3] We find that the num- 
bers of both merged and ejected particles increase be- 
yond At(e > 0.5)/t scc = 1, while most survived par- 
ticles have At(e > 0.5)/r soc < 1. A peak at 10 -5 for 
survived test particles represents all the particles with 
At{e > 0.5)/r sec < 10 -4 5 , and their forced eccentricities 
never exceed 0.5. By comparing each run with a figure 
similar to Fig. [3j we find that the following formula pre- 
dicts the fates of test particles well. 

At(e > 0.5)/r sec > 1 : merged, or ejected (5) 
< 1 : survived (6) 

3.3. Fates of Test Particles depending on Dynamical 
Outcomes of Jupiters 

In the previous sections, we have discussed the trends 
seen from all the simulations. We have shown that the 
fates of test particles are determined by secular inter- 
actions as well as close encounters with giant planets. 
Therefore, it is natural to expect that the fates of test 
particles depend on dynamical outcomes of giant planets 
as well. We explore this possibility here. 

By the end of the simulations, all forty systems have 
lost at least one Jupiter-mass planet, and thirteen sys- 
tems have lost two Jupiters. When we categorize the 
systems by the first-lost giant planet, 17/40, 15/40, 
and 8/40 systems have Jupiter- Jupiter collisions, Jupiter 
ejections, and Jupiter-Star mergers, respectively. As 
shown in Fig. for Jupiter-Jupiter collision cases, semi- 
major axes of survived Jupiter-like planets tend to be 
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tag (o/AU) Eccentricity 

Fig. 7. — Left: normalized semi-major axis distributions for the 
survived Jupiter-like planets. From top to bottom, Jupiter-Jupiter 
collision, Jupiter ejection, and Jupiter-Star merger cases. Right: 
normalized, cumulative distributions of corresponding eccentrici- 
ties. 
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Fig. 8. — Dependence of test particles' fates on dynamical out- 
comes of giant planets. Each panel shows the fates of test particles 
as a function of the initial semimajor axes of test particles normal- 
ized by the minimum semimajor axes of Jupiter-mass planets. Top, 
middle, and bottom panels correspond to the systems with colli- 
sions of two giant planets, ejections of giant planets, and mergers of 
giant planets with the central stars. It is apparent that most sur- 
vived test particles reside in Jupiter-Jupiter collision cases, while 
almost all test particles are gone for Jupiter-Star merger cases. 

similar to the initial values, and more than 90 % of these 
planets have eccentricities less than 0.5. On the other 
hand, for Jupiter ejection and Jupiter-Star merger cases, 
final semi-major axes of giant planets are more widely 
distributed, and less than 40 % of planets have eccentric- 
ities below 0.5. This plot indicates that, as expected, 
Jupiter- Jupiter collision cases are relatively dynamically 
quiet, while ejection and merger cases are more violent. 

In Fig. [BJ we show histograms similar to the bottom 
panel of Fig. U depending on Jupiters' dynamical out- 
comes. As shown in Fig. [7J when giant planets collide 
with each other, neither semi-major axes nor eccentric- 
ities of these planets change largely. Such evolution of 
giant planets should have the least significant effects on 
orbital stability of test particles. Indeed, the top panel 
shows that most "stable" test particles are in Jupiter- 
Jupiter collision cases. 

On the other hand, when a Jupiter merges with the 
central star (bottom panel), the planet sweeps through 
the area of test particles and removes most of them. 
There are several survived test particles as well. How- 



ever, as we can guess from their random distribution in 
atp.in/aj.min, ah" of them have high eccentricities. Note 
that, in this Jupiter-Star merger case, most test parti- 
cles also merge with the central stars rather than being 
ejected out of the systems. Thus, although giant planets 
approach test particles, the eccentricities of these parti- 
cles are excited not by close encounters, but by secular 
interactions. 

The middle panel shows the results of the Jupiter- 
ejection cases. Here, the pericenters of ejected Jupiters 
do not necessarily come close to the regions of test par- 
ticles. Therefore, some test particles can survive. How- 
ever, the majority of test particles either merge with the 
central stars or are ejected. 

For Jupiter- Jupiter collision, Jupiter ejection, and 
Jupiter-Star merger cases, the survival percentage of test 
particles are 61, 14.5, and 5.7%, respectively. Thus, our 
results suggest that, unless the evolution of giant plan- 
ets is either stable or quiet (e.g., Jupiter- Jupiter colli- 
sion cases), most test particles interior to Jupiters' orbits 
would be removed from the systems, even if their orbital 
radii are several times smaller than the pericenter dis- 
tance of the innermost survived Jupiters. Since systems 
which underwent dynamically violent evolutions tend to 
have giant pl anets on eccentric orbits, our results confirm 
the study bv lVeras fc Armitaed (|2006l ). and suggest that 
systems with eccentric Jupiters at a few AU may have 
lower chances of hosting terrestrial planets. We discuss 
this issue further in Section EOl and IH 

3.4. Effects of Dynamical Instability on Stable Regions 

Previous studies estimated the stability of hypotheti- 
cal te rrestrial planets in the HZs of the observed systems 
(e.g IGehman et al.lfl996l: I Jones et al.ll200ll: [Noble et all 
1200a IMenou fc Tabachnikl 120031: lAsghari et all I2004D . 
However, since these studies did not take account of the 
past evolution of the observed systems, they are likely 
to have overestimated the dynamically stable region of 
terrestrial planets. In the cases of Jupiter ejections and 
and Jupiter-Star mergers, Jupiters with the highest ec- 
centricities must have disappeared from the systems and 
these planets should have strongly destabilized the orbits 
of terrestrial planets through secular perturbations that 
were highly enhanced by the largest eccentricities. In 
this subsection, we show how the estimated stable region 
would change by ignoring the past orbital evolution. To 
study this, we add the same eleven test particles we used 
as initial conditions to the final state of giant planets for 
each simulation, and run the simulations for lOMyr. As 
expected, we find that many more test particles survive 
until the end of the simulations. 

Fig. compares the outcomes of two sets of simula- 
tions. The top panel is similar to the bottom panel of 
Fig. [H but shows the fates of test particles in terms of 
the initial semi-major axes of test particles normalized 
by the innermost Jupiter's semimajor axes at the end 
of the simulations at p ,i n /ajji n , instead of at p ,in/ ^J.min- 
The corresponding figure for the new set of simulations 
is shown in the bottom panel. 

In the top panel, we see the tail of ejected test particles 
with small a tp .in/ 'ajjin- This tail is a result of different 
normalization factors (i.e., ajji n instead of aj. m ; n ), and 
does not mean that there are many ejected particles inte- 
rior to the stable test particles. Most of the ejected test 
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Fig. 9. — The effects of dynamical instability of giant planets 
on the orbital stability of test particles. The top panel is similar 
to the bottom panel of Fig. [4] but the initial semimajor axis of 
a test particle is normalized by the final semimajor axis of the 
innermost giant planet. The bottom panel shows the corresponding 
results for simulations of 11 test particles added to the final states 
of giant planets. Many more test particles survive in the latter case 
compared to the former, which indicates that dynamical instability 
among giant planets indeed affects the stable region of the test 
particles. 

particles in this tail belong to the Jupiter-Star merger 
cases, where the survived giant planets often have large 
final semimajor axes (see Fig. EJ- 

While the overall shapes of these two distributions are 
similar, there is a significant increase in the fraction of 
stable test particles in the bottom panel, which confirms 
our expectation that past dynamical instabilities among 
giants strongly influence the orbital stability of test par- 
ticles. In other words, the stable regions of hypothetical 
terrestrial planets that are determined from the current 
distributions of giant planets would be largely overesti- 
mated. 

For the new simulations, the fractions of the survived 
test particles have increased for all kinds of giant planet 
evolutions, compared to the original simulatios. As ex- 
pected, the increase is most significant for Jupiter-Star 
merger cases, and very modest for Jupiter- Jupiter colli- 
sion cases. The large increase in the number of survived 
test particles for Jupiter-Star merger and Jupiter ejection 
cases indicates that, even when the stability analysis al- 
lows the existence of stable hypothetical terrestrial plan- 
ets, few such planets would exist interior to moderately 
to highly eccentric giant planets. However, we should 
note that our simulations don't take account of tidal cir- 
cularization effects, which may have saved some of these 
terrestrial planets. We briefly discuss this possibility in 
Section |U 

3.5. Comparisons with Observations 

In this subsection, we compare our numerical re- 
sults with the distribution of the observed exoplane- 
tary systems. Out of over 600 confirmed planetary 
systems, more than 100 are multiple-planet systems 
(http : //exoplane t . eu7| . We choose systems with both 
.Neptune- and Jupiter-mass planets, where we define 
these planets to have less and more masses than O.lMj, 
respectively. There are 13 systems like that including our 
solar system, and the semimajor axis ratios for Neptune- 
and Jupiter-mass planets are plotted in the middle panel 
of Fig. [TO] The eccentricities of giant planets in these 
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Fig. 10. — Comparison of our results shown in the top panel 
of Fig. [9] with the observed multiple-planet systems. For the mid- 
dle panel, we select all systems with both Neptune- and Jupiter- 
like planets, where these planets are defined to be less and more 
massive than O.lMj. The ratios of these planets' semimajor axes 
to the innermost Jupiter-like planets' semimajor axes are plotted. 
The observed systems fall in the stable test-particle area that is 
estimated from our simulations. For the bottom panel, Kepler's 
multiple-planet-candidate systems are shown. All systems with 
planets smaller and larger than 6Re are selected. The cyan circles 
show the smaller planets that have the radii comparable to Earth 
of < 2Re, while the pink circles represent the rest. Kepler planets 
tend to be much closely separated from each other, compared to 
the expectations from our numerical simulations. 

systems are lower than the other giant planets, and more 
than 90% of them have eccentricities less than 0.4. The 
preference for low eccentricities of giant planets in mixed- 
size multiple-planet systems agrees with the expectation 
from our simulations (see Section [3.3[) . 

In Fig. [TOj the circles on the same dotted line belong 
to the same system. We find that the confirmed sys- 
tems with both Neptune- and Jupiter-like planets are 
distributed within the range of survived test particles 
from our simulations (top panel). There are two circles 
exterior to this region, and beyond ajy/aj = 1. However, 
since our simulations don't have test particles outside the 
Jupiters' orbits initially, we don't intend to estimate the 
stable region beyond aisf/aj = 1 . 

We perform a similar comparison for Kepler's planet 
candidates as well, by choosing systems with Neptune- 
and Jupiter-size planets (smaller and larger than 6Re, 
respectively) . The limit is adopted so that it is consistent 
with the classification of Kepler's planetary candidates 
(e.g.. iBatalha et al.ll2012| ). Out of about 360 multiple- 
planet-candidate systems, 22 have both of these planets. 
As shown in the bottom panel of Fig.[TUl Kepler's planet 
candidates turn out to be more closely separated from 
each other than our simulations predict. However, the 
agreement is better for systems with smallest planets of 
< 2Re (cyan circles). There are a few potential reasons 
why the distributions of Kepler's planet candidates are 
less consistent with our simulations. First, it is possi- 
ble that most multiple-planet-candidate systems did not 
experience any dynamical instability. If this is the case, 
we would expect that these planetary candidates have 
small eccentricities. The future observations would help 
clarifying this possibility. Secondly, the agreement may 
improve if we take masses of small planets into account. 
We discuss this point further in the next section. Thirdly, 
our simulations are done with equal-mass giant planets, 
and thus are expected to place the lower limit on the 
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dynamically stable region of test particles. It is possible 
that many Kepler planet candidates are actually in the 
"stable" region. 

4. DISCUSSION AND CONCLUSIONS 

The recent observations discovered many multiple- 
planet systems, including the ones with multiple sub- 
Neptune-size planets on short-period orbits. Further- 
more, there are many giant-planet systems which allow 
the presence of terrestrial planets in HZs. However, it is 
still unclear whether there are many planetary systems 
in which giant and terrestrial planets coexist like the So- 
lar system. In fact, Kepler and HARPS surveys suggest 
that close-in sub-Neptune-size planets are generally not 
accompanied by a gas giant planet in the same system — 
~ 5 — 10 % of multiple-planet systems have such mixed- 
size planets. 

One of the potential causes that prevents the coex- 
istence of giant and terrestrial planets would be early 
evolution of giant-planet systems. The migration and 
dynamical instability of giant planets are two potentially 
importa nt effects on format i on an d survival of terrestrial 
planets. iVeras fe~Ar mitagc (2005) pointed out that when 
dynamical instability of giant planets occurs, terrestrial 
planets in the HZ undergo close encounters with the gi- 
ants and are ejected out of the system. 

In this paper, we have studied how dynamical insta- 
bility among gas giant planets affects the orbital stabil- 
ity of sub-Neptune-like planets interior to their orbits, 
particularly that of close-in terrestrial planets. We have 
performed N-body simulations of systems with three 
Jupiter-mass planets and test particles corresponding to 
terrestrial planets over 0.1 — 3 AU. The important results 
we have found are as follows: 

• Dynamical instability often leads to either an ejec- 
tion of a giant planet or a merger of such a planet 
with the central star. In these cases, test particles 
are removed over a broad range of semimajor axis 
a>tp,in down to several times smaller than the min- 
imum pericenter distance of the giants (7j,min 

(see 

Fig. [2} . This means that stable regions evaluated 
by the current distribution of giants is significantly 
overestimated. 

• A major cause of removal of test particles is sec- 
ular perturbations from giant planets. Eccentric- 
ity excitement of giants during their orbit crossing 
highly enhances the effect of secular perturbations. 
Since secular perturbations continuously increase 
eccentricities of the test particles without changing 
their semimajor axes, most of the removed particles 
merge with the host star rather than being ejected 
out of the system. 

• Even for orbit crossing cases of Cbtp,in/<lJmin > lj 
there are many mergers (top panel of Fig. 2]). Most 
of these test particles do not have actual close en- 
counters with giant planets, due to large mutual 
inclinations between the orbits, and/or tendency 
away from the pericenter misalignment (see Fig.[S])- 
Due to these effects, the merged and ejected par- 
ticles are divided around atp,in/ a 'J,min = E rather 
than a t p,in/qj,mm = 1- 



• Since most test particles are removed as a result 
of secular interactions but not of close scattering, 
we can semi-analytically estimate the dynamically 
unstable region of test particles (see Section 13.21) . 
We have found that the ratio of the total duration 
of the forced eccentricity ef orco d > 0.5 and the sec- 
ular time At(e > 0.5)/r scc is a good measure for 
this region (see Fig. . This formula allows us to 
estimate the unstable region of test particles just 
from the simulations of giant planets. 

When the orbital instability of gas giants does not oc- 
cur or only Jupiter- Jupiter collisions occur, secular per- 
turbations do not remove terrestrial planets over a wide 
range of orbital radii because the eccentricities of giants 
stay low. On the other hand, in the systems with giants 
of relatively high eccentricities, the secular perturbations 
are highly enhanced so that most terrestrial planets in- 
terior to the orbits of giants are removed except for ones 
very close to the host stars. Our work further indicates 
that the eccentricities of survived particles are not high, 
since secular perturbations sensitively depend on the ra- 
tio of semimajor axes of a test particle and giant planets. 
Our theoretical prediction of the orbital stable region is 
consistent with the observational data (see Section l3~5|) . 

Note that there are some caveats and limitations in 
our analysis. In principle, the lowest-order secular theory 
used for the semi-analytical formula cannot be applied to 
the high eccentricity cases. However, we have found that 
it still provides a fairly good agreement with the numeri- 
cal simulations. The dynamically stable region estimated 
from equal-mass Jupiters cases corresponds to the lower 
limit (jVeras fe Armitagdl2006l ). Self-gravity among ter- 
restrial planets may also inhibit instability, because grav- 
itational interactions among the massive particles that 
are perturbed by giant planets tend to align neighbor- 
ing particles with one another an d make their eccentric- 
ity evolution synchron ized (e.g., Ilto fe T anikawa 120011 : 
iLevison fc Agn or 2003). However, such an effect can also 
lead to mergers of all the planets with the central star all 
at once. Furthermore, we have implicitly assumed that 
dynamical instabilities among giant planets occur after 
the gas disk dissipates, although dynamical instabilities 
could occur once the d isk masses become com parable to 
planetary masses (e.g.. iMatsumura et al.|[2009f ). The gas 
disks also provide additional precessions of planetary or- 
bits and the dissipation of gas disks can move th e loca- 
tions of secular resonances (|Nagasawa et al.ll200Q[ ). Such 
an effect can influence the overall results, but is beyond 
the scope of this paper. The other effects which precess 
the planetary orbits and change the resonance locations 
include general relativistic corrections as well as rota- 
tional and tidal oblateness effects of a star and a planet. 
However, these effects are important only near the cen- 
tral stars, and thus are unlikely to change our conclusions 
drastically. 

Our results have some theoretical predictions and ob- 
servational implications as well: 

• Formation of very close-in sub-Neptunes by scatter- 
ing. We have not taken account of the tidal inter- 
actions between planets and the central star. With 
tidal effects, the highly eccentric giant planets can 
become close-in plan ets due to orbital decay and 
circularization (e.g. iFabrvckv fc Tremaind 120071 : 
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Fig. 11. — Comparison of various timescales. Blue solid 
line shows the secular timescale between an Earth-like planet at 
0.1 AU and a Jupiter-like planet at 1AU. Solid and dashed green 
curves show tidal circularization timescales for an Earth-like and 
a Jupiter-like planet at 0.1 AU by assuming tidal quality factors of 
100 and 10 5 , respectively. Solid and dashed red curves are corre- 
sponding orbital decay timescales for an Earth-like and a Jupiter- 
like planet, respectively. Secular timescale becomes less than cir- 
cularization timescale for an Earth-like planet for an eccetricity 
e > 0.9. From the orbital decay timescale, we find that such a 
circularized Earth-like planet could survive for up to 0.1 — lGyr. 



Wu et all [20071: iNagasawa et all [20081 : iNaoz et all 
2011tlWu fe Lithwickl 120111 ). Such a giant planet 



may remove closer-in terrestrial planets which sur- 
vived in our simulations. Alternatively, highly 
eccentric sub-Neptune-mass planets could also 
achieve very short period, circular orbits, if or- 
bital circularization timescales are short compared 
to secular timescales. We show comparisons of 
some relevant timescales in Fig. [TTJ For an Earth- 
like planet at 0.1 AU which is accompanied by a 
Jupiter-like planet at 1 AU, the secular timescale 
T sec .E-j becomes shorter than the tidal circular- 
ization timescale r Et E for an eccentricity e > 0.9. 
The corresponding orbital decay timescale r a , e for 
such a case can be up to 0.1 — 1 Gyr, which indicates 
that such a planet could survive for a typical age 
of a planet-hosting star. Similar tidal dissipation 
timescales for a Jupiter- like planet at 0.1 AU (not 
1 AU) are also plotted in dashed curves for com- 
parison. Note that the circularization timescale 
for a Jupiter-like planet becomes comparable to 
an Earth-like planet for tidal quality factors of 
Q = 10 5 and 100, respectively. 

Metallicity enhancement of the stellar surface due 
to the mergers. Since our study suggests that many 
terrestrial planets and/or planetary embryos could 
merge with the central star under violent evolu- 
tion of giant planets, it is informative to estimate 



whether such mergers would change the stellar 
metallicity significantly. For typical solar-like stars, 
a convective envelope accounts for a few percent of 
the total mass (about 0.02 M Q - 4 x 10 28 kg), and 
about 2 % of this mass is expected to be composed 
of metals (~ 8 x 10 26 kg). When we write the stel- 
lar metallicity as [Fe/H] = \og[Z/X] -log[Z Q /X Q ], 
and assume that the mergers of terrestrial planets 
only change the mass fraction of metals Z, the to- 
tal mass of ~ 35 Me is required to change the stel- 
lar metallicity by 0.1 from the solar value. Thus, 
a significant change of a stellar metallicity might 
be possible if there were multiple mergers of Super 
Earths. However, it could be difficult to tell the 
primordial metallicity from such a raised value by 
mergers. 

• Lack of sub-Neptunes around metal-rich stars. The 
giant planets have been found around solar-type 
stars with metallicities ranging from [Fc/H] ^ — 0.7 
to ~ 0.6, with preferences for metal-rich stars. 
On the other hand, sub-Neptune-like planets are 
not found around meta l-rich stars with [Fc/H] > 
0.2 (|Mavor et all 120111 ). Since higher metallici- 
ties imply higher formation rates of giant plan- 
ets, it is possible that sub-Neptune-mass plan- 
ets in metal-rich systems tend to be exposed to 
dynamical instabilities of giant planets. Then, 
dynamical instabilities would remove these sub- 
Neptune-mass planets as described in this paper, 
and they could be rarer among metal-rich sys- 
tems. Alternatively, a rec ent Monte Carlo popu- 
lation synthesis study by [Mordasi ni et al.l (|2012[ ) 
showed that their single planet formation model 
produced more sub-Neptune-mass planets around 
lower-metallicity stars. Currently, the metallic- 
ity distribution of Kepler's planet- hosting stars is 
poorly known, though there is an indication that 
these stars are metal rich compared to st ars with 
no planets (jSchlaufman fe Laughlinl 120111) . If the 
present metallicity trend were confirmed by fu- 
ture observations including Kepler, HARPS, and 
ESPRESSO, we would need to understand whether 
such a trend is due to natural planet formation 
or evolution involving dynamically unstable giant 
planets. 
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